home *** CD-ROM | disk | FTP | other *** search
/ SGI Developer Toolbox 6.1 / SGI Developer Toolbox 6.1 - Disc 4.iso / src / exampleCode / 3Dmodeling / curve.c < prev    next >
Encoding:
C/C++ Source or Header  |  1994-08-02  |  14.5 KB  |  528 lines
  1. /*
  2.  * Copyright 1992, 1993, 1994, Silicon Graphics, Inc.
  3.  * All Rights Reserved.
  4.  *
  5.  * This is UNPUBLISHED PROPRIETARY SOURCE CODE of Silicon Graphics, Inc.;
  6.  * the contents of this file may not be disclosed to third parties, copied or
  7.  * duplicated in any form, in whole or in part, without the prior written
  8.  * permission of Silicon Graphics, Inc.
  9.  *
  10.  * RESTRICTED RIGHTS LEGEND:
  11.  * Use, duplication or disclosure by the Government is subject to restrictions
  12.  * as set forth in subdivision (c)(1)(ii) of the Rights in Technical Data
  13.  * and Computer Software clause at DFARS 252.227-7013, and/or in similar or
  14.  * successor clauses in the FAR, DOD or NASA FAR Supplement. Unpublished -
  15.  * rights reserved under the Copyright Laws of the United States.
  16.  */
  17.  
  18. /* Tom Davis -- 1992 */
  19.  
  20. #include "3d.h"
  21. #include <stdio.h>
  22. #include <math.h>
  23. #define PI    3.1415926535
  24.  
  25. curve_t *newcurve()
  26. {
  27.     curve_t *c = (curve_t *)malloc(sizeof(curve_t));
  28.     c->type = 0;
  29.     c->data = c->cache = 0;
  30.     return c;
  31. }
  32.  
  33. typedef struct {
  34.     float (*x)(float), (*y)(float), (*z)(float);
  35.     float (*x1)(float), (*y1)(float), (*z1)(float);
  36.     float (*x2)(float), (*y2)(float), (*z2)(float);
  37. } analyticdata_t;
  38.  
  39. typedef struct {
  40.     float a[3], b[3], c[3], d[3];
  41. } cubicdata_t;
  42.  
  43. void freecurve(curve_t *c)
  44. {
  45.     (*c->fctns->cleanup)(c);
  46.     free(c);
  47. }
  48.  
  49. static void analyticpoint(curve_t *c, float t, float p[3])
  50. {
  51.     analyticdata_t    *cd = (analyticdata_t *)c->data;
  52.  
  53.     p[0] = (*cd->x)(t);
  54.     p[1] = (*cd->y)(t);
  55.     p[2] = (*cd->z)(t);
  56. }
  57.  
  58. static void analytictangent(curve_t *c, float t, float p[3])
  59. {
  60.     analyticdata_t    *cd = (analyticdata_t *)c->data;
  61.  
  62.     p[0] = (*cd->x1)(t);
  63.     p[1] = (*cd->y1)(t);
  64.     p[2] = (*cd->z1)(t);
  65. }
  66.  
  67. static void analyticaccel(curve_t *c, float t, float p[3])
  68. {
  69.     analyticdata_t    *cd = (analyticdata_t *)c->data;
  70.  
  71.     p[0] = (*cd->x2)(t);
  72.     p[1] = (*cd->y2)(t);
  73.     p[2] = (*cd->z2)(t);
  74. }
  75.  
  76. static void analyticcleanup(curve_t *c)
  77. {
  78.     free(c->data);
  79.     if (c->cache) free(c->cache);
  80. }
  81.  
  82. cblock_t analyticblock = {
  83.     analyticpoint,
  84.     analytictangent,
  85.     analyticaccel,
  86.     analyticcleanup
  87. };
  88.  
  89. curve_t *newanalyticcurve(float (*x)(float), float (*y)(float), float (*z)(float),
  90.         float (*x1)(float), float (*y1)(float), float (*z1)(float),
  91.         float (*x2)(float), float (*y2)(float), float (*z2)(float))
  92. {
  93.     curve_t    *cur = newcurve();
  94.     analyticdata_t    *cdata = (analyticdata_t *)malloc(sizeof(analyticdata_t));
  95.  
  96.     cur->type = CURVE_ANALYTIC;
  97.     cur->fctns = &analyticblock;
  98.     cur->data = (long *)cdata;
  99.     cdata->x = x; cdata->y = y; cdata->z = z;
  100.     cdata->x1 = x1; cdata->y1 = y1; cdata->z1 = z1;
  101.     cdata->x2 = x2; cdata->y2 = y2; cdata->z2 = z2;
  102.     return cur;
  103. }
  104.  
  105. static void cubicpoint(curve_t *c, float t, float p[3])
  106. {
  107.     cubicdata_t    *cd = (cubicdata_t *)c->data;
  108.     long i;
  109.  
  110.     for (i = 0; i < 3; i++)
  111.     p[i] = cd->d[i] + t*(cd->c[i] + t*(cd->b[i] + t*cd->a[i]));
  112. }
  113.  
  114. static void cubictangent(curve_t *c, float t, float p[3])
  115. {
  116.     cubicdata_t    *cd = (cubicdata_t *)c->data;
  117.     long i;
  118.  
  119.     for (i = 0; i < 3; i++)
  120.     p[i] = cd->c[i] + t*(2.0*cd->b[i] + t*3.0*cd->a[i]);
  121. }
  122.  
  123. static void cubicaccel(curve_t *c, float t, float p[3])
  124. {
  125.     cubicdata_t    *cd = (cubicdata_t *)c->data;
  126.     long i;
  127.  
  128.     for (i = 0; i < 3; i++)
  129.     p[i] = 2.0*cd->b[i] + 6.0*t*cd->a[i];
  130. }
  131.  
  132. static void cubiccleanup(curve_t *c)
  133. {
  134.     free(c->data);
  135.     if (c->cache) free(c->cache);
  136. }
  137.  
  138. cblock_t cubicblock = {
  139.     cubicpoint,
  140.     cubictangent,
  141.     cubicaccel,
  142.     cubiccleanup
  143. };
  144.  
  145. curve_t *newcubiccurve(float a[3], float b[3], float c[3], float d[3])
  146. {
  147.     curve_t    *cur = newcurve();
  148.     cubicdata_t    *cdata = (cubicdata_t *)malloc(sizeof(cubicdata_t));
  149.  
  150.     cur->type = CURVE_CUBIC;
  151.     cur->fctns = &cubicblock;
  152.     cur->data = (long *)cdata;
  153.     copy3(a, cdata->a);
  154.     copy3(b, cdata->b);
  155.     copy3(c, cdata->c);
  156.     copy3(d, cdata->d);
  157.     return cur;
  158. }
  159.  
  160. /* For the Bezier, B Spline, and Cardinal splines, the following
  161.  * code converts them to cubic splines so that they can be expanded
  162.  * as polynomials.  The conversion is simply a change of basis
  163.  * accomplished by a matrix multiplication.  See, for example, 
  164.  * "Geometric Modeling", by Michael Mortensen.
  165.  */
  166.  
  167. static float bez[4][4] = {
  168.     -1.0, 3.0, -3.0, 1.0,
  169.     3.0, -6.0, 3.0, 0.0,
  170.     -3.0, 3.0, 0.0, 0.0,
  171.     1.0, 0.0, 0.0, 0.0
  172. };
  173.  
  174. curve_t *newbeziercurve(float p[3], float q[3], float r[3], float s[3])
  175. {
  176.     curve_t    *cur = newcurve();
  177.     cubicdata_t    *cdata = (cubicdata_t *)malloc(sizeof(cubicdata_t));
  178.     long    i;
  179.  
  180.     cur->type = CURVE_CUBIC;
  181.     cur->fctns = &cubicblock;
  182.     cur->data = (long *)cdata;
  183.     for (i = 0; i < 3; i++) {
  184.     cdata->a[i] = bez[0][0]*p[i]+bez[0][1]*q[i]+bez[0][2]*r[i]+bez[0][3]*s[i];
  185.     cdata->b[i] = bez[1][0]*p[i]+bez[1][1]*q[i]+bez[1][2]*r[i]+bez[1][3]*s[i];
  186.     cdata->c[i] = bez[2][0]*p[i]+bez[2][1]*q[i]+bez[2][2]*r[i]+bez[2][3]*s[i];
  187.     cdata->d[i] = bez[3][0]*p[i]+bez[3][1]*q[i]+bez[3][2]*r[i]+bez[3][3]*s[i];
  188.     }
  189.     return cur;
  190. }
  191.  
  192. static float bsp[4][4] = {
  193.     {-1.0/6.0, 3.0/6.0, -3.0/6.0, 1.0/6.0},
  194.     {3.0/6.0, -6.0/6.0, 3.0/6.0, 0.0},
  195.     {-3.0/6.0, 0.0,  3.0/6.0,  0.0},
  196.     {1.0/6.0, 4.0/6.0, 1.0/6.0, 0.0}
  197. };
  198.  
  199. curve_t *newbsplinecurve(float p[3], float q[3], float r[3], float s[3])
  200. {
  201.     curve_t    *cur = newcurve();
  202.     cubicdata_t    *cdata = (cubicdata_t *)malloc(sizeof(cubicdata_t));
  203.     long    i;
  204.  
  205.     cur->type = CURVE_CUBIC;
  206.     cur->fctns = &cubicblock;
  207.     cur->data = (long *)cdata;
  208.     for (i = 0; i < 3; i++) {
  209.     cdata->a[i] = bsp[0][0]*p[i]+bsp[0][1]*q[i]+bsp[0][2]*r[i]+bsp[0][3]*s[i];
  210.     cdata->b[i] = bsp[1][0]*p[i]+bsp[1][1]*q[i]+bsp[1][2]*r[i]+bsp[1][3]*s[i];
  211.     cdata->c[i] = bsp[2][0]*p[i]+bsp[2][1]*q[i]+bsp[2][2]*r[i]+bsp[2][3]*s[i];
  212.     cdata->d[i] = bsp[3][0]*p[i]+bsp[3][1]*q[i]+bsp[3][2]*r[i]+bsp[3][3]*s[i];
  213.     }
  214.     return cur;
  215. }
  216.  
  217. static float card[4][4];
  218.  
  219. static void makecardinal(float a)
  220. {
  221.     card[0][0] = card[1][3] = card[2][0] = -a;
  222.     card[0][3] = card[2][2] = a;
  223.     card[2][1] = card[2][3] = card[3][0] =
  224.     card[3][2] = card[3][3] = 0.0;
  225.     card[3][1] = 1.0;
  226.     card[0][1] = 2.0 - a;
  227.     card[0][2] = -2.0 + a;
  228.     card[1][0] = 2.0 * a;
  229.     card[1][1] = -3.0 + a;
  230.     card[1][2] = 3.0 - 2.0 * a;
  231. }
  232.  
  233. curve_t *newcardinalcurve(float p[3], float q[3], float r[3], float s[3])
  234. {
  235.     curve_t    *cur = newcurve();
  236.     cubicdata_t    *cdata = (cubicdata_t *)malloc(sizeof(cubicdata_t));
  237.     long    i;
  238.     static long inited = 0;
  239.  
  240.     if (inited == 0) {
  241.     inited = 1;
  242.     makecardinal(.5);
  243.     }
  244.  
  245.     cur->type = CURVE_CUBIC;
  246.     cur->fctns = &cubicblock;
  247.     cur->data = (long *)cdata;
  248.     for (i = 0; i < 3; i++) {
  249.     cdata->a[i] =card[0][0]*p[i]+card[0][1]*q[i]+card[0][2]*r[i]+card[0][3]*s[i];
  250.     cdata->b[i] =card[1][0]*p[i]+card[1][1]*q[i]+card[1][2]*r[i]+card[1][3]*s[i];
  251.     cdata->c[i] =card[2][0]*p[i]+card[2][1]*q[i]+card[2][2]*r[i]+card[2][3]*s[i];
  252.     cdata->d[i] =card[3][0]*p[i]+card[3][1]*q[i]+card[3][2]*r[i]+card[3][3]*s[i];
  253.     }
  254.     return cur;
  255. }
  256.  
  257. #define SSTEPS 400
  258.  
  259. /* curvelength() breaks the curve into SSTEPS steps in parameter
  260.  * space, and then approximates it with a broken line segment
  261.  * whose length is calculated.  It would be better to do this
  262.  * adaptively, with more steps where the curvature is higher.
  263.  */
  264.  
  265. float curvelength(curve_t *c)
  266. {
  267.     float    len = 0.0, *cache;
  268.     float    p[3], q[3];
  269.     long    i;
  270.  
  271.     (*c->fctns->point)(c, 0.0, p);
  272.     cache = (float *)malloc((1+SSTEPS)*sizeof(float));
  273.     cache[0] = 0.0;
  274.     c->cache = (long *)cache;
  275.     for (i = 1; i <= SSTEPS; i++) {
  276.     (*c->fctns->point)(c, (float)i/(float)(SSTEPS), q);
  277.     len += dist3(p, q);
  278.     cache[i] = len;
  279.     copy3(q, p);
  280.     }
  281.     return len;
  282. }
  283.  
  284. /* curvestep returns in p the point on the curve n steps into it --
  285.  * i.e. at the point (float)(n/nsegs).
  286.  */
  287.  
  288. static float curvestep(curve_t *c, long n, long nsegs, float len, float p[3])
  289. {
  290.     long    i;
  291.     float    limit, l1, l2;
  292.     float    p1[3], q1[3], val;
  293.     float    *cache = (float *)c->cache;
  294.     static long nlast = -1, ilast;
  295.     static float    plast[3];
  296.     static float    vallast;
  297.  
  298.     if (n == nlast) {
  299.     p[0] = plast[0]; p[1] = plast[1]; p[2] = plast[2];
  300.     return vallast;
  301.     }
  302.     if (n == 0) {
  303.     nlast = -1;
  304.     (*c->fctns->point)(c, 0.0, p);
  305.     return 0.0;
  306.     }
  307.     if (n == nsegs) {
  308.     nlast = -1;
  309.     (*c->fctns->point)(c, 1.0, p);
  310.     return 1.0;
  311.     }
  312.     limit = len*(float)n/(float)nsegs;
  313.     for (i = 1; i < SSTEPS; i++)
  314.     if (cache[i] > limit) {
  315.         val = ((float)(i-1) +
  316.         (limit-cache[i-1])/(cache[i]-cache[i-1]))/((float)SSTEPS);
  317.         (*c->fctns->point)(c, val, p);
  318.         nlast = n;
  319.         plast[0] = p[0]; plast[1] = p[1]; plast[2] = p[2];
  320.         vallast = val;
  321.         return val;
  322.     }
  323. }
  324.  
  325. pwlin_t    *pwlinfromcurve(curve_t *cur, long n)
  326. {
  327.     pwlin_t    *pwl = newpwlin();
  328.     float    *f, len;
  329.     long    i;
  330.  
  331.     pwl->n = n+1;
  332.     pwl->verts = f = (float *)malloc((n+1)*3*sizeof(float));
  333.     len = curvelength(cur);
  334.     for (i = 0; i <= n; i++)
  335.     (void) curvestep(cur, i, n, len, &f[i*3]);
  336.     return pwl;
  337. }
  338.  
  339. void makesmoothpwlworm(curve_t *cur, pwlin_t *pw, long nsegs, void (*savefunc)())
  340. {
  341.     makepwlworm(cur, pw, nsegs, smoothsave);
  342.     doverts();
  343.     doquads(savefunc);
  344. }
  345.  
  346. /* Makes a "worm" by moving a Fresnel frame along the curve.  The
  347.  * tangent points along the curve, the normal is the derivative of
  348.  * the tangent, and the binormal is perpendicular to both.  See any
  349.  * elementary book on differential geometry for more information.
  350.  * The code in the following routine, makeworm(), is similar.
  351.  */
  352.  
  353. void makepwlworm(curve_t *cur, pwlin_t *pw, long nsegs, void (*savefunc)())
  354. {
  355.     static float lasttan[3] = { 1, 0, 0, };
  356.     static float lastnorm[3] = { 0, 1, 0, };
  357.  
  358.     float    slen;
  359.     float    p[3], t, tan[3], norm[3], binorm[3];
  360.     float    tan1[3], norm1[3], binorm1[3];
  361.     float    p0[3], p1[3], p2[3], p3[3], pp[3];
  362.     float    n0[3];
  363.     long    i, j;
  364.     float    s, c, tt, s1, c1, t1, dot;
  365.     long    nsides = pw->n - 1;
  366.  
  367.     slen = curvelength(cur);
  368.     for (i = 0; i < nsegs; i++) {
  369.     t = curvestep(cur, i, nsegs, slen, p);
  370.     (*cur->fctns->tangent)(cur, t, tan);
  371.     if (tan[0] == 0.0 && tan[1] == 0.0 && tan[2] == 0.0) {
  372.         error("makepwlworm: zero tangent");
  373.         copy3( lasttan, tan );
  374.     }
  375.     copy3( tan, lasttan );
  376.     (*cur->fctns->accel)(cur, t, norm);
  377.     if (norm[0] == 0.0 && norm[1] == 0.0 && norm[2] == 0.0) {
  378.         error("makepwlworm: zero acceleration");
  379.         copy3( lastnorm, norm );
  380.     }
  381.     copy3( norm, lastnorm );
  382.     crossprod(tan, norm, norm);
  383.     crossprod(tan, norm, binorm);
  384.     normalize(tan);
  385.     normalize(norm);
  386.     normalize(binorm);
  387.     t = curvestep(cur, i+1, nsegs, slen, pp);
  388.     (*cur->fctns->tangent)(cur, t, tan1);
  389.     if (tan1[0] == 0.0 && tan1[1] == 0.0 && tan1[2] == 0.0) {
  390.         error("makepwlworm: zero tangent");
  391.         copy3( lasttan, tan1 );
  392.     }
  393.     copy3( tan1, lasttan );
  394.     (*cur->fctns->accel)(cur, t, norm1);
  395.     if (norm1[0] == 0.0 && norm1[1] == 0.0 && norm1[2] == 0.0) {
  396.         error("makepwlworm: zero acceleration");
  397.         copy3( lastnorm, norm1 );
  398.     }
  399.     copy3( norm1, lastnorm );
  400.     crossprod(tan1, norm1, norm1);
  401.     crossprod(tan1, norm1, binorm1);
  402.     normalize(tan1);
  403.     normalize(norm1);
  404.     normalize(binorm1);
  405.     for (j = 0; j < nsides; j++) {
  406.         tt = pw->verts[3*j + 2];
  407.         s = pw->verts[3*j + 1];
  408.         c = pw->verts[3*j];
  409.         t1 = pw->verts[3*j + 5];
  410.         s1 = pw->verts[3*j + 4];
  411.         c1 = pw->verts[3*j + 3];
  412.         p0[0] = c*norm[0] + s*binorm[0] + tt*tan[0] + p[0];
  413.         p0[1] = c*norm[1] + s*binorm[1] + tt*tan[1] + p[1];
  414.         p0[2] = c*norm[2] + s*binorm[2] + tt*tan[2] + p[2];
  415.         p1[0] = c1*norm[0] + s1*binorm[0] + t1*tan[0] + p[0];
  416.         p1[1] = c1*norm[1] + s1*binorm[1] + t1*tan[1] + p[1];
  417.         p1[2] = c1*norm[2] + s1*binorm[2] + t1*tan[2] + p[2];
  418.         p2[0] = c1*norm1[0] + s1*binorm1[0] + t1*tan1[0] + pp[0];
  419.         p2[1] = c1*norm1[1] + s1*binorm1[1] + t1*tan1[1] + pp[1];
  420.         p2[2] = c1*norm1[2] + s1*binorm1[2] + t1*tan1[2] + pp[2];
  421.         p3[0] = c*norm1[0] + s*binorm1[0] + tt*tan1[0] + pp[0];
  422.         p3[1] = c*norm1[1] + s*binorm1[1] + tt*tan1[1] + pp[1];
  423.         p3[2] = c*norm1[2] + s*binorm1[2] + tt*tan1[2] + pp[2];
  424.  
  425.         if (samepoint(p0, p1) || samepoint(p1, p2))
  426.         perpnorm(p0, p2, p3, n0);
  427.         else
  428.         perpnorm(p0, p1, p2, n0);
  429.  
  430.         scalarmult(-1.0, n0, n0);
  431.         m_xformpt(p0, p0, n0, n0);
  432.         m_xformptonly(p1, p1);
  433.         m_xformptonly(p2, p2);
  434.         m_xformptonly(p3, p3);
  435.  
  436.         (*savefunc)(ADD_QUAD, n0, p0, n0, p1, n0, p2, n0, p3);
  437.     }
  438.     }
  439. }
  440.  
  441. void makeworm(curve_t *cur, float radius, long nsides, long nsegs, void (*savefunc)())
  442. {
  443.     float slen;
  444.     float p[3], t, tan[3], norm[3], binorm[3];
  445.     float tan1[3], norm1[3], binorm1[3];
  446.     float p0[3], p1[3], p2[3], p3[3], pp[3];
  447.     float n0[3], n1[3], n2[3], n3[3];
  448.     long i, j;
  449.     float s, c, s1, c1, dot;
  450.  
  451.     slen = curvelength(cur);
  452.     for (i = 0; i < nsegs; i++) {
  453.     t = curvestep(cur, i, nsegs, slen, p);
  454.     (*cur->fctns->tangent)(cur, t, tan);
  455.     if (tan[0] == 0.0 && tan[1] == 0.0 && tan[2] == 0.0) {
  456.         error("makeworm: zero tangent");
  457.         return;
  458.     }
  459.     (*cur->fctns->accel)(cur, t, norm);
  460.     if (norm[0] == 0.0 && norm[1] == 0.0 && norm[2] == 0.0) {
  461.         error("makeworm: zero acceleration");
  462.         return;
  463.     }
  464.     crossprod(tan, norm, norm);
  465.     crossprod(tan, norm, binorm);
  466.     normalize(norm);
  467.     normalize(binorm);
  468.     t = curvestep(cur, i+1, nsegs, slen, pp);
  469.     (*cur->fctns->tangent)(cur, t, tan1);
  470.     if (tan1[0] == 0.0 && tan1[1] == 0.0 && tan1[2] == 0.0) {
  471.         error("makeworm: zero tangent");
  472.         return;
  473.     }
  474.     (*cur->fctns->accel)(cur, t, norm1);
  475.     if (norm1[0] == 0.0 && norm1[1] == 0.0 && norm1[2] == 0.0) {
  476.         error("makeworm: zero acceleration");
  477.         return;
  478.     }
  479.     crossprod(tan1, norm1, norm1);
  480.     crossprod(tan1, norm1, binorm1);
  481.     normalize(norm1);
  482.     normalize(binorm1);
  483.     for (j = 0; j < nsides; j++) {
  484.         s = sin(2*j*PI/nsides);
  485.         c = cos(2*j*PI/nsides);
  486.         if (j+1 == nsides) {
  487.         s1 = sin(0.0);
  488.         c1 = cos(0.0);
  489.         } else {
  490.         s1 = sin(2*(j+1)*PI/nsides);
  491.         c1 = cos(2*(j+1)*PI/nsides);
  492.         }
  493.         n0[0] = c*norm[0] + s*binorm[0];
  494.         n0[1] = c*norm[1] + s*binorm[1];
  495.         n0[2] = c*norm[2] + s*binorm[2];
  496.         p0[0] = radius*n0[0] + p[0];
  497.         p0[1] = radius*n0[1] + p[1];
  498.         p0[2] = radius*n0[2] + p[2];
  499.         n1[0] = c1*norm[0] + s1*binorm[0];
  500.         n1[1] = c1*norm[1] + s1*binorm[1];
  501.         n1[2] = c1*norm[2] + s1*binorm[2];
  502.         p1[0] = radius*n1[0] + p[0];
  503.         p1[1] = radius*n1[1] + p[1];
  504.         p1[2] = radius*n1[2] + p[2];
  505.         n2[0] = c1*norm1[0] + s1*binorm1[0];
  506.         n2[1] = c1*norm1[1] + s1*binorm1[1];
  507.         n2[2] = c1*norm1[2] + s1*binorm1[2];
  508.         p2[0] = radius*n2[0] + pp[0];
  509.         p2[1] = radius*n2[1] + pp[1];
  510.         p2[2] = radius*n2[2] + pp[2];
  511.         n3[0] = c*norm1[0] + s*binorm1[0];
  512.         n3[1] = c*norm1[1] + s*binorm1[1];
  513.         n3[2] = c*norm1[2] + s*binorm1[2];
  514.         p3[0] = radius*n3[0] + pp[0];
  515.         p3[1] = radius*n3[1] + pp[1];
  516.         p3[2] = radius*n3[2] + pp[2];
  517.  
  518.         m_xformpt(p0, p0, n0, n0);
  519.         m_xformpt(p1, p1, n1, n1);
  520.         m_xformpt(p2, p2, n2, n2);
  521.         m_xformpt(p3, p3, n3, n3);
  522.  
  523.         (*savefunc)(ADD_QUAD, n0, p0, n1, p1, n2, p2, n3, p3);
  524.     }
  525.     }
  526. }
  527.  
  528.